Influence of external flows on crystal growth: numerical investigation 
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We use a combined phase-field/lattice-Boltzmann scheme [D. Medvedev, K. Kassner, Phys. Rev. 
E 72, 056703 (2005)] to simulate non-facetted crystal growth from an undercooled melt in external 
flows. Selected growth parameters are determined numerically. 

For growth patterns at moderate to high undercooling and relatively large anisotropy, the values of 
the tip radius and selection parameter plotted as a function of the Peclet number fall approximately 
on single curves. Hence, it may be argued that a parallel flow changes the selected tip radius and 
growth velocity solely by modifying (increasing) the Peclet number. This has interesting implications 
for the availability of current selection theories as predictors of growth characteristics under flow. 

At smaller anisotropy, a modification of the morphology diagram in the plane undercooling versus 
anisotropy is observed. The transition line from dendrites to doublons is shifted in favour of dendritic 
patterns, which become faster than doublons as the flow speed is increased, thus rendering the basin 
of attraction of dendritic structures larger. 

For small anisotropy and Prandtl number, we find oscillations of the tip velocity in the presence 
of flow. On increasing the fluid viscosity or decreasing the flow velocity, we observe a reduction in 
the amplitude of these oscillations. 

PACS numbers: 47.54.-r, 47.20.Hw, 02.70.-c, 68.70.+W. 
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I. INTRODUCTION 

Crystal growth from the melt or from solution almost 
never occurs in convection-free conditions. Notwith- 
standing this fact, models of solidification often focus, 
when dealing with aspects of morphological stability and 
pattern formation, on situations where convection is ab- 
sent pi |E IE • Rather than by the negligibility of con- 
vective effects, such a choice is generally motivated by the 
difficulty of including them and the argument that the 
basic prototypes of patterns appear and may be studied 
without convection. This in turn has led to the tailor- 
ing of experiments on dendritic growth, in which it was 
explicitly tried to avoid disturbing flows 0, IE IE • 

When convection was taken into account in calcula- 
tions, it was usually in simplified geometries or in con- 
ditions enabling simplification of the model [E IE fioj : 
often the moving-boundary aspect of the problem was 
neglected ^E El Attempts to model the full prob- 
lem including convection and the motion of the liquid- 
solid interface were essentially made only in cases where 
the deflection of the interface remained relatively small 
0,0] . The state of the art a decade ago may be summa- 
rized roughly by saying that pattern formation in crystal 
growth could be well simulated either for the solid, treat- 
ing the free-boundary problem in its full complexity, or 
for the liquid, obtaining the convection roll pattern with 
good accuracy by use of simplifying assumptions for in- 
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terface motion. 

With the advent of efficient phase-field techniques , 
the solution of the moving-boundary problem became 
simpler, and first simulations of convection in dendritic 
growth were performed in diverse geometries and with 
both imposed flows and natural convection |l6L IrH Hg| . 
In these, the Navier-Stokes equations were solved by 
standard numerical approaches implementing the partial 
differential equations either within a finite-element or a 
finite-difference scheme. It was then a natural idea to 
supplement the efficient approach to interface motion by 
an efficient non-iterative method for flow simulation, the 
lattice-Boltzmann scheme. This approach was pioneered 
by Miller et al. [3 HE EH, and a slightly different vari- 
ant, the advantages of which will be discussed in Sec.llTl 
was developed by ourselves EE Ei| . 

Numerical studies of pattern formation solving the 
combined free-boundary and flow problems will be useful 
in guiding the development of analytic selection theory 
for dendritic growth and other growth modes in the pres- 
ence of convection. At present, there is a well-developed 
theory for purely diffusion-limited dendritic growth, both 
in two and three dimensions [H EE EE E3 EE E& EE 
l3l| . It provides an analytic demonstration of the exis- 
tence of a discrete set of needle crystal solutions to the 
model equations and shows that the fastest of these so- 
lutions is linearly stable. 

Acceptance of this microscopic solvability theory has 
not been uncontroversial, as there is no clear agreement 
of its predictions referring to crystalline anisotropy with 
experimental results ^ij^ . On the other hand, the mathe- 
matical statements of the theory can hardly be disputed. 
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Therefore, the existence of a needle crystal solution is a 
fact and its stability shows that it is an attractor of the 
dynamics. In principle, its basin of attraction might be 
so small as to render it irrelevant experimentally. But 
this is essentially excluded by numerical simulations that 
have shown both in two [s^lls^ l and three dimensions 
that the dynamical state of the system more or less auto- 
matically approaches the prediction of selection theory. 
In two dimensions at least, this is also true for the den- 
drite decorated with side branches. Numerics and theory 
agree with each other so that experimental discrepancies 
are most likely due to the fact that anisotropics in real 
crystals do not correspond to the model expressions used 
in the theory or else that additional effects interfere which 
are absent or controlled in the simulations (for example, 
thermal fluctuations not considered in the deterministic 
theory might affect the operating point of the dendrite 
at low undercoolings [3(|). 

One of the more surprising predictions of microscopic 
solvability theory was the nonexistence of dendrites in 
the absence of any kind of anisotropy, be it that of sur- 
face tension or of interface kinetics. Due to the nature 
of the theoretical approach, which is singular perturba- 
tion theory about an Ivantsov parabola or paraboloid 
[33, such a statement can hold only for needle crystals 
with a shape close to the (exact) solution of which they 
are supposed to be small perturbations. And it turned 
out later that indeed steady-state crystal growth at zero 
anisotropy is possible, but only with a shape that is far 
from an Ivantsov solution. These new structures were 
called doublons [34] in two dimensions and selection the- 
ory has been developed for them as well [3^. Since they 
continue to exist at finite anisotropy, there is a coexis- 
tence range with dendrites, which means that there are 
two attractors of the dynamics. The standard argument 
is then that the faster of the two morphologies will win, 
which in the analytically tractable case is the doublon, 
whenever it exists. 

Large scale two-dimensional structures consist of ar- 
rays of dendrites or of doublons evolving in a noisy envi- 
ronment via side branching and/or tip splitting processes. 
A theoretical description of the resulting dendrite and 
seaweed morphologies, based on scaling arguments and 
selection theory [39j gives a kinetic phase diagram in the 
parameter space of undercooling versus surface tension 
anisotropy. All of the analytic predictions mentioned so 
far refer to diffusion-limited growth only. 

A first attempt to extend selection theory to situa- 
tions with a flow was made by Bouissou and Pelce 01 
and there were a number of less rigorous approaches to 
the problem as well (references are given in [23j). How- 
ever, while one experiment seems to support this theory 
plf. another one contradicts it |4^. Moreover, the the- 
ory is based on a linearization of the basic equations, an 
approach that has been found not to always be reliable 
|43|| . Clearly, more numerical or experimental data are 
needed to both check the existing theories and guide fur- 
ther theoretical development. The purpose of this article 



is to provide first elements of these data usin g ou r com- 
bined phase-field lattice-Boltzmann approach |22l |23| . A 
successful selection theory on the microscale (constructed 
on the basis of these and similar data) would then yield 
useful information for more applied work on macrosegre- 
gation and related questions. 

The paper is organized as follows. In Sec. 2, we discuss 
the basic model equations and describe the method of 
their numerical implementation. In Sec. 3, we consider 
the influence a parallel flow on the selection of growth 
velocity and tip radius, whereas in Sec. 4, changes in 
the position of the morphology transition from dendrites 
to doublons induced by the flow are discussed. Some 
concluding remarks are given in Sec. 5. 

II. MODEL EQUATIONS AND 
IMPLEMENTATION 

For simplicity, we consider a symmetric model with 
equal thermal diffusivities of the solid and liquid phases, 
expecting it to display all the qualitative features of the 
more general case. Moreover, since we wish to confront 
our simulations with theoretical results on different mor- 
phologies, we restrict ourselves to two-dimensional sys- 
tems here. So far, there are very few results on non- 
dendritic structures in three dimensions. 

The well-accepted sharp-interface description of non- 
facetted crystal growth from a supercooled melt in the 
presence of a fluid flow with velocity U then starts from 
the following set of bulk and interface equations: 

dtU + UVu = DW 2 u, 

n V = Dn-(V«|.-Vtt|,), (1) 
u t = -d{9)/R K -f3n- V. 

Herein, u = c p (T — T m )/L with c p denoting the heat 
capacity and L the latent heat, both per unit volume, is 
a nondimensionalized temperature, T being the temper- 
ature at the considered position and T m the bulk melt- 
ing temperature. D is the thermal diffusivity, n the 
local normal to the liquid-solid interface (pointing from 
the solid into the liquid) and V the interface velocity. 
The subscripts of the temperature gradients in the sec- 
ond equation refer to the solid and the liquid sides of 
the interface, respectively. d(9) is the capillary length, 
related to the orientation-dependent surface tension "f(9) 
by d{6) = (do/70) [7(f) + l"{9)}. 9 is the angle between 
the interface normal and some fixed direction (usually 
identified with the x axis of the coordinate system), 70 
the average of the interface tension over all orientations, 
and do = joT m c p / L 2 the similarly averaged capillary 
length. Rk is the local radius of curvature, (3 the ki- 
netic coefficient. In principle, (3 is an orientation depen- 
dent quantity just as the surface tension; however, we 
will restrict ourselves to the case of fast attachment ki- 
netics here, meaning that (3 becomes negligible. This 
implies certain constraints on the choice of parameters of 
the phase-field model (see below). 
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The boundary condition for the normal velocity is usu- 
ally referred to as Stefan condition, while the last equa- 
tion in describing the interface temperature, is the 
Gibbs-Thomson relation (for j3 = 0) with kinetic correc- 
tion (for (3 ^ 0). 

At infinity, the temperature in the solid approaches 
T m , corresponding to w = 0, whereas in the liquid, it 
takes on some value < T m , corresponding to u = 
—A. The quantity A is denoted as the nondimcnsional 
undercooling. 

In order not to complicate things unnecessarily, we as- 
sume the melt to be an incompressible Newtonian fluid, 
governed by the appropriate version of the Navier-Stokes 
equations, supplemented by boundary conditions at the 
interface 



d t U + UVU 
V -u 



VP 



0, 
0. 



(2) 



where equal mass densities p have been assumed in the 
two phases, v is the kinematic viscosity, and P denotes 
the pressure of the liquid. Ui is the flow velocity at the 
interface. The boundary conditions correspond to no-slip 
conditions for the tangential velocity and describe the 
fact that for equal densities of the two phases, the liquid 
is neither sucked towards nor rejected from the interface; 
hence, its normal velocity at the interface is zero in the 
laboratory frame (the rest frame of the solid). 

It is convenient to use dimensionless variables in the 
analysis of the growth process. The tip radius can be 
nondimensionalized using the capillary length as R = 
R/d , whereas flow and growth velocities become nondi- 
mensional via normalization with the characteristic ve- 
locity given by the ratio of the thermal diffusion coef- 
ficient and the average capillary length, that is U = 
Ud /D and V = Vd /D. 

The material properties are characterized by the 
anisotropy of the surface energy or rather that of the 
capillary length and by the Prandtl number Pr = v j ' D\ 
the flow is characterized by the Reynolds number Re = 
UR/v. 

We employ a combined phase- field/lattice-Boltzmann 
scheme where solidification is simulated with the phase- 
field model of Karma and Rappel [35l |44| , and the flow 
of the liquid as well as convective and diffusive heat 
transport are modeled using a lattice-Boltzmann (LB) 
method. This means that the actual equations simulated 
are not those given above but a phase-field approxima- 
tion to the interface dynamics (involving a finite-width 
transition region between the solid and its melt) and a set 
of kinetic equations that are asymptotically equivalent to 
the Navier-Stokes and advection-diffusion equations. 

In particular, the phase-field equations read 

T(6)8til> = (^-Au(l-^ 2 )) (1-V 2 ) 
+ V - (W 2 (6)VtP) 



- d x {w{9)w\e)d y ^) 

+ d y (W{6)W'(6)d x i>) , 



(3) 



ftw + UVu = DV 2 u+ -d t h(ip). 



The value i/j = 1 of the phase-field variable is chosen to 
represent the solid phase, whereas ip = — 1 corresponds to 
the liquid phase. W(9) is an anisotropic interface width, 
t(9) a relaxation time, and 9 = a,rcta,n(dy4> / d x ip) is the 
angle between the normal on a level set of ip and the 
x axis. For the level set given by rp = 0, this angle is 
the same as the angle 6 in Eq. otherwise it provides 
an extension of the definition of the latter into the bulk. 
h(ip) describes the coupling of the diffusion field to the 
phase field via latent heat production. This function was 
chosen as h(tp) = -ijj, which appears to be computation- 
ally most efficient |35| - 

Via a suitable asymptotic expansion, the equations of 
the sharp-interface model can be derived [35j, with 
the following expressions for the capillary length and ki- 
netic coefficient: 



d(6) 
0(6) 
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\JW{9) 



d 2 e W{9)) , 

K + JF W 2 {9) 
21 t(9) 
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These equations, first given by Karma and Rappel |44j . 
have been shown to be equivalent to a second-order ac- 
curate standard asymptotic approximation [3^, 0] . 
Requiring 



W = W A(9), 



t A 2 (9), A = 



2IDt q 



(K + JF)W 2 



we can impose a vanishing kinetic coefficient [44). 
For h(ijj) = ip, the values of the phase-field spe- 
cific coefficients are / = 2V2/3, J = 16/15, K = 
\/2(Hf - £§ln2), and F = v^hi2 M, Ei We 
use the anisotropy function 



leading to 



A(6)= 'A = 1 + — cos 40. 
7o 15 



<io(l — a cos A9) , 



(4) 



which is the usual model expression exhibiting four-fold 
symmetry. 

Moreover, we set to = 1, Wo = 1. 

The equation for the phase-field tp was discretized on 
a uniform spatial lattice with grid spacing Ax = 0.4, 
and it was solved using the explicit Euler method with 
constant time step At in the range 0.008-0.016. 

Both the flow of the liquid and the heat transport are 
simulated using the LBGK method (see ^s main 

variables are one-particle distribution functions /j. de- 
fined on the nodes of a regular spatial lattice. Differ- 
ent labels k correspond to different velocities from a 
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fixed finite set. In the two-dimensional model used here, 
we employ nine velocities, namely, c = (0,0), = 
(cos((fc - l)7r/2),sin((fc - l)-ir/2))8x/5t for k = 1...4, 
and c fc = V2(cos((fc - 1/2)tt/2), sin((fc - 1/2)-k/2))8x/ 8t 
for k = 5 ... 8. Here, 5x is the grid spacing, equal for 
both directions, 5t is the time step. The effect of making 
the velocities proportional to 5x/St is that nonzero ve- 
locities lead to nearest neighbour and next-nearest neigh- 
bour sites of the square lattice in one time-step, i.e., only 
lattice point positions appear in the dynamics, no inter- 
polations are necessary. 

Inside the LBE part of our simulations, the grid spac- 
ing and time step are both formally rescaled to one, which 
is the reason why we have used a different notation for 
them here from that in the phase-field part of the sim- 
ulation (5x and St vs. Ax and At), although they are 
actually the same "material" quantities. 

The evolution equation for f k is 



f k (t + St,x + c k St)=f k (t,x) + 



f eq 
Jk 



fk 



St. (5) 



Distribution functions are advected (first term on the 
r.h.s.) and undergo a relaxation to equilibrium val- 
ues f^ q which are, as usual, taken to be expansions of 
Maxwellians up to second order in the fluid velocity U 



c fe -U , (c fc -U) 2 



2c! 



El 

2c? 



(6) 



with c s having the physical meaning of an isothermal 
sound velocity. The local fluid density is given by p — 
J2 fk = Yl fk 9 > tne now velocity is U = J2 /fcCfc/p, and 

the weight coefficients are wq = 4/9, iUi_4 = 1/9, u>5_s = 
1/36. This form of the equilibrium distribution functions 
ensures mass and momentum conservation and provides 
the correct form of the momentum flux tensor [4g, [42| • 

Performing a Chapman-Enskog expansion, one can de- 
rive from (jSJ) the continuity and Navier-Stokes equations 
[46| . with kinematic viscosity v = c 2 (jf — St/2). The 
isothermal sound velocity is c s = Sx/V3St, for small 
flow velocities the fluid is almost incompressible (effects 
of compressibility are proportional to U 2 /c 2 ). 

The influence of the growing pattern on the fluid flow 
was simulated as proposed in |l7l Hsj . An additional 
dissipative force was introduced in partially filled regions 

where 5=2.757 and <j> — (1 + ip)/2 is the solid fraction. 
This provides the correct velocity boundary conditions at 
the diffuse interface (see i.e., the sharp-interface 

limit of the velocity boundary conditions of Eq. (J2J is 
correctly reproduced. The value of the constant g was 
obtained in |l7l ITsI ] by an asymptotic analysis of plane 
flow past the diffusive interface. 

The action of forces on a liquid was simulated by 
the exact difference method of 48] . The term Af k = 



fl q {p,\5 + AU) - /fc*(/o,U) is added to the r.h.s. of 
eq. (J5J, where AU = FSt/p is the velocity change due 
to action of force F during time step St. This form of 
the change of distribution functions is exact in the case 
where the distribution is equilibrium before the action of 
the force (then after the action the distribution remains 
equilibrium), whence the name of the method. In the 
case of a non-equilibrium initial state, this method is ac- 
curate to second order in AU. It is simple enough and 
valid for arbitrary lattices and any dimension of space. 

The heat transport equation in was treated in a 
very similar way via the introduction of nine additional 
distribution functions N k (t, x). An extensive discussion 
of these algorithmic details is presented in |2^| . 

In order to give a general impression of the type of re- 
sults obtainable with these simulations, we display Figs.^ 
and|21 showing the steady states of a dendritic and a dou- 
blon pattern, respectively. 



FIG. 1: Symmetric needle crystal, i.e., dendritic pattern in 
antiparallel flow. Growth parameters: A = 0.75, a = 0.45, 
U — 0.01. The capillary length is do = 0.185, the measured 
growth velocity V = 0.0451, leading to a diffusion length of 
8.2 (20.5 lattice units). The flow pattern is indicated by the 
streaks outside of the crystal. Numerical grid size: 700 x 
1400 corresponding to 1513.5 do X 3027.0 do. 

A few remarks may be in order comparing our model 
with its main predecessor as given by Miller et al. 
[l9l I20L l2lj. There are obvious similarities, but our ap- 
proach is simpler in the lattice-Boltzmann part and has 
better convergence in the phase-field part. The model 
discussed in l20l 1211 ] is four-dimensional and uses 24 
velocities. We have a two-dimensional model with 2 x 
9 velocities and our collision matrix is simpler. So the 
lattice-Boltzmann part of our model is faster in two di- 
mensions, which is the only case considered here. 

Concerning the phase- field part, our approach includes 
the Karma model which has been shown to be quantita- 
tive at much larger interface thicknesses than preceding 
alternatives. The phase-field model used by Miller et al. 
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FIG. 2: Asymmetric needle crystal, (half of a) doublon pat- 
tern. Same growth parameters as in Fig. except that 
a = 0.3. The capillary length and system size are the same as 
in Fig.^as well, the measured growth velocity is V — 0.0402, 
leading to a diffusion length of 9.2, or 23 lattice units. 

has not been demonstrated to have any of the advan- 
tages of the thin-interface asymptotics. Its quantitative 
accuracy might be challenged on the basis of the same 
objections as that of the original Kobayashi model |4j|. 

Therefore, with the same accuracy prescribed, we ex- 
pect the phase-field part of our model to be much more ef- 
ficient (because it converges to the correct sharp-interface 
limit at much smaller system sizes) than that of the 
model used in [T^.l20ll2l| and the lattice-Boltzmann part 
to be slightly more efficient. 

III. SELECTION OF GROWTH PARAMETERS 

The growth of a single needle crystal in parallel flow 
was simulated for fixed surface tension anisotropy and a 
range of undercoolings and flow velocities. To investigate 
the effects of parallel flow on the growth characteristics, 
needle crystals grown without flow were used as initial 
configurations. Reaching the steady state in the absence 
of flow took between 90000 and 300000 time steps. 

After loading the values of the temperature and phase 
fields, the flow was initialized with boundary condi- 
tions of constant flow velocity perpendicular to the up- 
per boundary and zero velocity gradients at the lower 
boundary, while the side boundaries were made reflect- 
ing. First, the flow was allowed to evolve with a fixed 
configuration of the solid, and the relative velocity error 
was calculated at each time step as 

TT _ Y.\Ux-Ux\ + \Uy-Uy\ 

E\U X \ + \Uy\ ■ 

Here, U refers to the flow velocity at the current, U to 



that at the preceding time step, the summation is over all 
grid nodes. The convergence condition was U err < 10~ 5 . 

Then the growth of the pattern was "switched on" and 
continued until a steady, i.e., constant-velocity, state was 
reached. In the range of undercoolings 0.72 < A < 0.76, 
this took on the order of 150000 time steps. 

The numerical grid in these runs had a size of 700 
x 1400, corresponding to between 505 g?o x 1010 do and 
2014 d x 4028 d . For the smallest measured velocities, 
the diffusion length remained below 250 lattice units, for 
the largest one, it was about 15 lattice units. Therefore, 
in all cases the system size was large enough to consider 
finite-size effects negligible, in particular in view of the 
fact that the computational domain corresponded to half 
the system size only (see Figs. ^ and |5J). 

All the simulations discussed in this section were done 
either until convergence of the pattern to a steady state 
was achieved or such a steady state was found to be 
unattainable - below we report on the appearance of 
oscillatory states in certain parameter regions. Only 
then growth characteristics such as the growth veloc- 
ity were measured, i.e., care was taken to avoid tran- 
sient states providing only information about temporary 
growth characteristics. 

Computed values of the reduced tip radius R and se- 
lection parameter a = 2Ddo/R 2 V = 2/R 2 V are plotted 
versus the growth Peclet number Pe = RV/2D = RV/2 
in Fig. and Fig. \Qp for dendrites (single symmetric 
fingers). In the figure, the anisotropy of surface stiffness 
is a — 0.75, the range of nondimensional initial under- 
coolings A = c p (T m — T 00 )/L extends from 0.4 to 0.8, 
and the reduced flow velocity U = Udo/D is typically 
chosen between and 0.32 (0, 0.01, 0.02, 0.04 etc.). One 
can see that for each of the two data sets most of the 
points fall onto a unique curve. Minor deviations appear 
mainly for small Prandtl numbers and large flow veloci- 
ties. The range of Reynolds numbers investigated in this 
data set was < Re < 5.6. It is possible to define a 
relative Reynolds number based on the flow speed in the 
reference system attached to the surface of the growing 
crystal. This Reynolds number, which was never zero, of 
course, extended up to « 7.1. 

To obtain the tip radius, we fitted a parabola to an 
extended region about the tip. This procedure does not 
yield an approximation to the inverse curvature at the 
tip itself, because due to anisotropy, the shape deviates 
from a parabola close to the tip [5(|. Rather it de- 
fines the tip radius by a global parabolic envelope as- 
sociated with energy conservation; the Peclet number 
Pe calculated from this tip radius corresponds to the 
Peclet number used in selection theory and in the ab- 
sence of flow reduces to the Ivantsov value, defined by 

A = V^Pe exp(Pe) erfc (VPej . 

The tip radius initially decreases with the Peclet num- 
ber but later begins to increase again (for this anisotropy, 
the minimum is at about Pe = 1), while the selec- 
tion parameter a decreases in the whole range of Peclet 
numbers investigated. That the tip radius increases for 
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and flow velocity at infinity, is in itself a nontrivial task. 
In limiting cases (small external flow speed), it may be 
approximated by the value obtained for an Ivantsov-type 
solution of an Oseen approximation to the problem with 
flow. 

According to selection theory for the purely diffusion- 
limited case, we should expect a to become independent 
of the Peclet number for small anisotropy and small un- 
dercooling. The latter condition can be relaxed - as 
long as Pea 1 / 2 <C 1, the standard result of selection the- 
ory V ~ (D/do) a 7 / 4 Pe 2 , continues to hold for a model 
anisotropy of the type (@J. However, due to computa- 
tional limitations, this limit is difficult to access, hence 
neither of these conditions is well satisfied in Fig. |3t>, 
where a = 0.75 and Pe varies through 1. The opposite 
limit of large Peclet number is also known analytically 
|5l|: the selection parameter should vary, for fixed small 
anisotropy, proportional to 1/Pe 2 . Moreover, it is possi- 
ble to evaluate the predictions of solvability theory jE^] 
numerically for arbitrary Peclet numbers. Formally, this 
can be done for arbitrary values of the anisotropy param- 
eter a - three examples are exhibited in Fig.Q]- but the 
theory should not be expected to yield good results for 
anisotropics that are not sufficiently small. 



FIG. 3: Dependence of the tip radius (a) and selection pa- 
rameter a (b) on the Peclet number for dendrites. Anisotropy 
parameter a — 0.75. Each symbol corresponds to the under- 
coolings and Prandtl numbers given in the legend and a value 
of the velocity U of the imposed flow in the range from to 
0.32. Larger values of U correspond to larger Peclet numbers. 
The dashed line in (b) is a fit to /(Pe) = a/ (1 + bPe) 2 , the 
dashed line in (a) is computed as g(Pe) = 1/ [/(Pe) Pe]. From 
the fit, a = 0.178, 6 = 0.841. 



large undercooling can be easily understood: in the limit 
A = 1 , the solution should approach a planar front with 
R = oo. This argument is made more quantitative below. 

From the theoretical point of view, the most interest- 
ing feature of these results is the existence of a master 
curve, on which the data fall for a wide range of pa- 
rameters. For this feature (if it holds generally) allows 
us to use the theory of dendritic growth without con- 
vection to make predictions of selected velocities and tip 
radii in the presence of a forced flow. In the absence 
of flow, the growth Peclet number depends only on the 
undercooling. As soon as flow is introduced, the Peclet 
number depends both on the undercooling and the ve- 
locity of the imposed flow. What Fig. |3Jd then tells us is 
essentially that no matter how we produce a given Peclet 
number, we should expect the same selected value of a at 
fixed anisotropy. Hence, if we change both the flow veloc- 
ity and the undercooling in a way that keeps the Peclet 
number constant, the growth speed and shape remain 
unchanged. This means that the case with flow can be 
mapped to the case without flow. Of course, the problem 
of determining the Peclet number, for given undercooling 
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FIG. 4: Behaviour of the selection parameter as a function 
of the Peclet number according to the linearized version of 
selection theory as given in |52| . The dash-dotted lines are 
fits with the same analytic expression as in Fig. [3Jj. 

Comparing Fig. [3] with Fig. 0] we see that the numer- 
ical results and the predictions from solvability theory 
show the same general trend of a decreasing with increas- 
ing Peclet number and that the order of magnitude of our 
a values is the same as for a = 0.75 in solvability theory. 
However, the selection parameter decreases much faster 
with increasing Peclet number than in the theory (note 
the different scales of the Pe axes). This may not be too 
surprising - after all the theory is made only for q« 1, 
and while it has been claimed to be quantitatively not 
too bad for a up to 0.5 or 0.6 this claim referred to 
the small Peclet number case. Since the Peclet number 
appears only at next-to-leading order in the small param- 
eter of the theory (-\/o : ) 1 it is perhaps not unreasonable 
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to assume that the dependence of selected characteristics 
on it is described less accurately within the theory than, 
say, the anisotropy dependence at Pe = 0. 

Because the theory predicts the limiting behaviours of 
the selection parameter at small and large Peclet num- 
bers, it is tempting to try to capture the behaviour at 
intermediate Pe by a simple interpolating function. The 
simplest rational function approaching a constant value 
for Pe — > at finite slope and being proportional to 1 /Pe 2 
at large Pe is /(Pe) = a/ (1 4- &Pe) . Fits with this func- 
tion work pretty well for both our numerical data and 
the results from selection theory as is demonstrated in 
Figs. 03 and 01 Because R = l/(<rPe), this also yields a 
description of the tip radii in Fig. |2t as well as a quan- 
titative estimate for the tip radius behavior at large Pe 
(~ Pe6 2 /a). The value of a gives the limit of the selec- 
tion parameter for Pe = 0, which in our case is about 
30% below the value obtained from selection theory (for 
an illegitimately large value of a). 

Note that given the graphs of a(Pe) and -R(Pe), we 
could obtain a similar representation for the growth ve- 
locity simply by plotting V = 2/(oR 2 ) or, even sim- 
pler, V = 2aPe 2 . Hence, the limiting growth velocity 
for Pe — ► oo is V = 2a/b 2 . 

Results for doublons, i.e., two asymmetric fingers with 
a liquid-filled channel between them growing together, 
are presented in Fig. [S] The surface stiffness anisotropy 
is 0.30 in this case; the undercooling ranges from 0.77 to 
0.85. For the reduced flow velocity the same range from 
to 0.32 was taken as for the symmetric finger, whereas 
the Reynolds numbers extended only up to 2.1, as the set 
of considered viscosities did not contain values as small 
as those of Fig. OH 

The "tip" radius was measured by fitting a parabola 
to the exterior shape of the doublon, that is, only points 
much farther from the central channel than the two tips 
of the pattern were used in the fitting procedure. Since 
this procedure depends also on the cutoff value defining 
which part of the shape is "exterior" and which one is 
"interior" , we do not expect a similar accuracy for this 
radius as in the case of dendritic patterns. Moreover, 
the total range of radii displayed is about a factor of six 
smaller than in Fig.yp, which contributes to making the 
results appear much noisier than those for the dendritic 
pattern. 

Nevertheless, while the characteristic length scale of 
the doublon may not display the same clear-cut univer- 
sality as that of the dendrite, a unique dependence of the 
selection parameter on the Peclet number is clearly visi- 
ble in Fig. 03. A fit to the same rational function /(Pe) 
as in the dendritic case is possible, but less trustworthy 
than for the dendrites, as the range of accessible Peclet 
numbers is smaller. Moreover, its theoretical foundation 
is weaker than for dendrites, as selection theory for dou- 
blons does not yet seem to have been worked out in the 
limit of large Peclet number. 

Finally, it should be mentioned that the introduction of 
an external flow may lead to a loss of stability of steady- 
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FIG. 5: Dependence of the tip radius (a) and selection pa- 
rameter a (b) on the Peclet number for doublons. Anisotropy 
parameter a — 0.3. Appearance of the same symbol several 
times means different values for U (in the same range as in 
Fig. 1) at the same pair of values for the undercooling and 
the Prandtl number. The dashed lines are obtained by fitting 
as in Fig. 01 which yields a = 0.0876, b = 0.699. 



state growth and result in instationary patterns. With 
small anisotropy and Prandtl number, oscillations of the 
tip velocity are observed. This observation may relate 
to the prediction by the selection theory presented in 
[40| that above a certain flow velocity no steady-state 
solutions will be possible any more. Increase of the fluid 
viscosity and/or decrease of flow velocity damps these 
oscillations as shown in Fig. H3 

From the existence of these oscillations, it may be 
concluded that there are parameter regions (attained 
for given anisotropy on decreasing the Prandtl number) 
where the simple picture discussed above breaks down. 
Selected growth characteristics will then not depend on 
the growth Peclet number and the anisotropy parameter 
alone. For these oscillations are not predicted by solv- 
ability theory without flow, hence the simple mapping to 
this theory is not feasible anymore, and an extension of 
selection theory such as the one given in [^j j but prefer- 
ably on a more rigorous basis, becomes necessary. 
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The rather strong flow provokes oscillations of the tip speed 
with large amplitude for small viscosity (Pr = 1.78) and small 
amplitude for large viscosity (Pr = 5.00). 



dritic and doublon morphologies |3J,[53 (at zero flow), it 
should be kept in mind that these older numerical results 
refer to the one-sided model whereas here we consider 
the symmetric model. As it turns out, the transition line 
is shifted to higher values of the anisotropy (e.g. a be- 
tween 0.20 and 0.25 at A = 0.7 instead of a sa 0.12), 
which seems plausible, because the added diffusion in 
the solid tends to reduce anisotropy-induced tempera- 
ture differences imposed at the interface. In fact, phase- 
field simulations of the symmetric model (in the absence 
of flow) by Tokunaga and Sakaguchi [54J also exhibit 
this shift to higher anisotropics. Their calculations were 
done in a channel that is narrow in comparison with our 
system width, so they introduced an intermediate mor- 
phology between doublons and dendrites, two competing 
Saffman- Taylor [55t l56[ like fingers (ST) [the existence 
range of which should vanish for infinite system size]. 
For A = 0.7, they find the transition from doublons to 
ST at a ss 0.24, that from ST to dendrites at a w 0.32, 
which agrees well with our result. 



IV. MORPHOLOGY DIAGRAM 

Previously, a kinetic phase diagram was obtained the- 
oretically |23 in the case of purely diffusive growth, dis- 
tinguishing four morphologies: compact dendritic struc- 
tures at large anisotropy and not too large undercooling, 
compact seaweed patterns at large undercooling (and not 
too large anisotropy), noise-dominated fractal dendritic 
and seaweed morphologies at sufficiently small anisotropy 
and undercooling, respectively. The transition lines be- 
tween the different morphologies and their nature (as 
first- or second-order kinetic phase transition or cross- 
over) were determined analytically under certain limit 
assumptions. Regarding the compact growth morpholo- 
gies, it was stated that doublons cease to exist for larger 
anisotropies, but when they exist, they are faster than 
dendrites. In principle, the latter exist at all nonzero 
anisotropies, but they are overtaken and thus overgrown 
by doublons in the region of coexistence (hence the tran- 
sition from compact dendrite to compact seaweed would 
be first order, because both morphologies coexist above 
a certain undercooling). 

How an imposed external flow may influence the differ- 
ent growth patterns is interesting and largely unexplored. 
We have already shown that doublons survive in a shear 
flow H3,|2^|, a somewhat counterintuitive result. 

In the present work, we investigate the morphology 
diagram for growth in a parallel flow imposing a num- 
ber of different flow velocities, with a particular view to 
the positions of the transition lines between doublon and 
dendrite growth. 

Figure gives an overview of the measured morphol- 
ogy diagram (actually a small section only of the entire 
plane undercooling versus anisotropy) for the purely dif- 
fusive case and two different flow velocities. In relating 
this to previously measured transitions between the den- 
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FIG. 7: Morphology diagram displaying the predominance 
of dendrites or doublons at different flow speeds. Triangles 
correspond to dendrites either being the only morphology or 
the faster one, inverted triangles correspond to doublons be- 
ing faster. There are three sizes of symbols. The smallest 
triangles refer to the case without flow described in current 
analytic theories. Medium-size triangles are for a flow speed 
U = 0.01; big ones for a flow speed U = 0.04. The gen- 
eral trend is that with increasing flow speed dendrites invade 
the original domain of doublons. There are three points at 
A = 0.67, where the simulation gave the same velocities for 
both structures to three significant digits. 

Because we simulate only one of the two fingers of 
a doublon, imposing mirror symmetry about the sys- 
tem boundary (see Fig. |2J), our calculation suppresses 
possible antisymmetric instabilities of a doublon, e.g., 
instabilities, where one finger gets ahead of the other. 
However, there is some evidence [34l that on increase 
of the anisotropy parameter doublons normally get un- 
stable by dynamical unbinding of the two fingers, which 
move apart and become independent dendrites. This un- 
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binding instability is symmetric and would not be missed 
by our approach. All our statements about existence of 
doublons are, of course, not affected by the possibility of 
an unstable mode not taken into account. And finally, 
we base our assertions about the predominance of one 
of the two growth modes on comparisons of the veloci- 
ties of both, which will come out correctly of the com- 
putation with the imposed symmetry. The worst that 
could happen is that a doublon found to be faster than 
a dendrite at the same parameter values is unstable with 
respect to an antisymmetric perturbation, in which case 
the dendritic morphology would survive, if it is stable. 
Such a scenario is not very likely, given the fact that our 
doublons, whenever they were faster than the associated 
dendrites, exhibited closely-spaced tips , corresponding to 
the predictions of selection theory |38( . 

The case of purely diffusive growth is depicted in Fig.0 
by the smallest symbols. Triangles with their tips point- 
ing upward correspond to dendrites, inverted triangles to 
doublons. On increase of the reduced flow velocity U 
to 0.01, denoted by larger triangles, dendrites become 
faster than doublons at several combinations of under- 
cooling and anisotropy. The largest triangles in Fig. [7| 
correspond to a velocity of U = 0.04. They demonstrate 
how the region where dendrites are faster than doublons 
increases with increasing flow velocity. 

It should be noted that according to the analytic the- 
ory [3^ for the purely diffusive case doublons would al- 
ways be faster than dendrites at coexistence. Dendrites 
would dominate only where doublons did not exist. This 
is not quite true at the large anisotropy values consid- 
ered here. For example, at the point a = 0.5, A = 0.67, 
where we have put a symbol denoting dendritic growth, 
doublons exist, too, but are slower than dendrites. Never- 
theless, it is rather remarkable that external flow can lead 
to dendrites becoming sufficiently fast to outrun doublons 
in an extended range of parameters. 

Note that the morphology diagram should actually be 
displayed in three dimensions, as it is spanned by the 
three variables a, A, and U. We circumvent the need for 
a genuine 3D representation by taking different symbol 
sizes to represent different flows, as only few flow velocity 
values could be studied. 

That the presence of a parallel fluid flow in general 
favours dendrites over doublons is an additional possible 
reason for the difficulty to obtain doublons in experi- 
ments, the main reason of course being that in experi- 
ments the value of the undercooling is usually so small 
that one is far from the existence region of doublons. 



Experimental approaches to produce doublons in crystal 
growth either had to use artifices to obtain effectively 
vanishing anisotropy in the growth plane l57| or led to 
the observation of transient doublons only |58Ll59l| - these 
were however true 3D structures. 

V. CONCLUSIONS 

In summary, we use a previously proposed combined 
phase- field/lattice-Boltzmann scheme to simulate den- 
dritic growth from a supercooled melt in external coun- 
terflows directed parallel to the growing needle crystal. 
Several regions of the morphology diagram in the space 
spanned by the anisotropy parameter, the nondimen- 
sional undercooling and the nondimensional flow velocity 
were explored. 

For dendrites at moderate to high undercooling and 
high anisotropy, we found that the values of tip radius 
and selection parameter, and hence of the growth veloc- 
ity, depend on the growth Peclet number only, not on 
the undercooling and flow velocity separately. Hence, 
it may be argued that the essential effect of a paral- 
lel flow, at least in a certain part of parameter space, 
is to change the selected tip radius and growth velocity 
solely by modifying (increasing) the Peclet number. In 
this region, selection theory for the purely diffusive case 
is applicable, the main task being to determine the re- 
lationship between undercooling, imposed flow velocity 
and the growth Peclet number. For doublons a similar 
dependence for the selection characteristics was obtained. 

Incorporation of genuine flow effects into selection the- 
ory does become necessary as the anisotropy and Prandtl 
number become small, when tip oscillations take over and 
the steady state either ceases to exist or becomes unsta- 
ble. Increase of the fluid viscosity and/or decrease of flow 
velocity is observed to damp down these oscillations. 

For smaller anisotropy, an interesting phenomenon is 
observed. The growth velocity for dendrites increases 
faster than for doublons with increase of the flow ve- 
locity (at the same undercooling and anisotropy). For 
some parameters, dendrites become faster, hence, exter- 
nal flow can lead to morphology transitions and change 
the kinetic phase diagram. 
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